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Streamwise and quasi-streamwise elongated structures have been shown to play a signif- 
icant role in turbulent shear flows. We model the mean behavior of fully turbulent plane 
Couette flow using a streamwise constant projection of the Navier Stokes equations. This 
results in a two-dimensional, three velocity component {2D/3C) model. We first use a 
steady state version of the model to demonstrate that its nonlinear coupling provides the 
mathematical mechanism that shapes the turbulent velocity profile. Simulations of the 
2D/3C model under small amplitude Gaussian forcing of the cross-stream components 
are compared to DNS data. The results indicate that a streamwise constant projection 
of the Navier Stokes equations captures salient features of fully turbulent plane Couette 
flow at low Reynolds numbers. A system theoretic approach is used to demonstrate the 
presence of large input-output amplification through the forced 2D/3C model. It is this 
amplification coupled with the appropriate nonlinearity that enables the 2D/3C model to 
generate turbulent behaviour under the small amplitude forcing employed in this study. 



1. Introduction 

The Navier Stokes (NS) equations provide a complete dynamical description of the 
three velocity components and pressure for simple canonical flows under the sole model- 
ing assumption that all important physical phenomena are captured by these equations. 
Unfortunately, these infinite dimensional algebraically constrained equations are ana- 
lytically intractable. They have however, been extensively studied computationally and 
numerical solutions do exist. For plane Couette flow, the first numerical solution was 



computed by Nagata (1990). Gibson et al. (2009) provide a detailed discussion of other 



work related to a full range of numerical plane Couette flow solutions. Ever increasing 
computing power will continue to allow progress toward understanding these local prop- 
erties. However, a full mathematical understanding of NS even in simple parallel flow 
configurations remains elusive, hence considerable effort has been applied to the search 
for more analytically-tractable flow models. 

In contrast, the Linearized Navier Stokes (LNS) equations can be analyzed using well 
developed tools from linear systems theory. For wall bounded shear flows, one particu- 
lar property of the LNS that has been extensively studied is disturbance amplification, 
e.g. iFarrelll (|1988|); iGustavssonl (|1991|); iReddy & Henningsonl (|1993|); iFarrell & loannou 



2 D. F. Gayme, B. J. McKeon, A. Papachristodoulou, B. Bamieh & J. C. Doyle 

(19936); Bamieh & Dahleh (2001); Jovanovic fc Bamieh| f2005). Large disturbance am- 
plification is common in these flows because the linear operators governing them are 
non-normal, (i.e., the operator, A, is such that A* A 7^ AA^). 

The LNS are thought to capture the energy production of the full nonlinear system. 



Henningson & Reddy (1994) showed that non- normality and linear mechanisms are nec- 



essary conditions for subcritical transition to turbulence and it is widely believed that 
energy ampliflcation is due to coupling terms that remain in linearized models, see for 
example ( jTrefethen et a/.||1993 ). In smooth wall-bounded shear flows the linear coupling 
between the Orr-Sommerfeld and Squire equations associated with nonzero spanwise 
wave number has been shown to be required for the generation of the wall layer streaks 
that are necessary to maintain turbulence (Butler fc FarreII||1993 Kim & Lim |2000} . In 
this context, the term "streak" describes the "well-deflned elongated region of spanwise 
alternating bands of low and high speed fluid" (Wale ffe et a/. [1991 ') . The LNS have also 
been used by Jovanovic fc Bamieh (2001) to predict certain second-order statistics of 
turbulent channel flow. The above results and a host of others illustrate the power of the 
LNS as a model for wall-bounded shear flows. There is however, one fundamental flow 
feature that linear models are unable to capture; the change in the mean velocity proflle 
as the flow transitions from laminar to turbulent. In addition, linear analysis can only 
give local information regarding the full (nonlinear) system. 

Empirical models have been shown to be useful in capturing key aspects of many 
flows. For example. Proper Orthogonal Decomposition (POD) has been successfully used 



to construct accurate low dimensional ordinary differential equation models, e.g. (Lumley 



1967 Smith et al. 2005). However the preceding analysis utilizes existing experimental 



or numerical data, a limitation also applicable to eddy viscosity models. In general, data- 
driven or heuristic models can be said to suffer from a lack of connection, of varying 
degree, to the governing equations of the problem. 

The model studied herein is an attempt to merge the beneflts of studying a physics- 
based set of equations, such as NS, with the analytical tractability of a simplifled model, 
such as the LNS. It is developed based on the assumption that certain aspects of fully 
developed turbulent flow can be reasonably modeled as homogeneous in the streamwise 
direction, here denoted "streamwise constant". The idea that a streamwise constant 
model is suflicient to capture mean proflle changes from laminar to turbulent is strongly 
supported by the work of|Reddy fc loannou (2000), who showed that nonlinear interaction 
between the ikx^kz) = (0, iA/") modes, where the /c's are the streamwise and spanwise 
wave numbers, is the primary factor in determining the turbulent mean velocity proflle 
in Couette flow. Further, as was discussed in Orlandi fc Jimenez ( 1994| ), this type of 
model may be adequate to capture many of the effects associated with the generation of 
turbulent wall fric tion. A 2^D mode l along similar lines has also been developed for the 



viscous wall layer; Tullis fc Pollard (1993) for example use such a model to study flow 



over riblets in this region. The physical and analytical basis for assuming homogeneity 
in the streamwise direction is discussed in the following section. 



1.1. Streamwise Coherence 



A growing body of work supports the notion that turbulence in wall-bounded shear flows 
is characterized by dynamically signiflcant coherent structures, particularly features with 
streamwise and quasi-streamwise alignment. Near- wall streaks (Kline et al. 1967), for 
example, have been shown to play a key role in energy production through the 'near-wall 
autonomous cycle' discussed by Waleffe (1990); Hamilton et al. (1995); Waleffe (1997); 



Jimenez fc Pinelli (1999). This cycle is generally agreed to be an important mechanism 
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in determining the low-order sta tistics of turbule nt flows in the buffer region and viscous 
sublayer, i.e. ^ 30 ( Schoppa fc Hussain||2QQ2" ). 

More recent high Reynolds number studies have focused on the identification and 
characterization of streamwise coherence in the core, i.e. Kim fc Adrian| ( 1999 ) ; |Morrison| 



et al. ( 2QQ4J ; ,Guala et a/. ( 2QQ6J ; Hut chins fc Marusic| (2QQ7a). These motions have been 



called large and very large scale motions (respectively LSM and VLSMs). They appear 
to have a similar signature to the near-wall streaks (Hut chins fc Marusic||2QQ76 Chung 
Ifc McKeon 2010), but tend to be longer in extent, from one to ten times the outer 



length scale, 5. There is experimental evidence to suggest that at high Reynolds numbers 
(for example Re^ > 7300), VLSMs contain more energy than the near- wall structures 
( [Morrison et a/.|[2004[ Hutchins fc Marusic||2007a[6 ). In turbulent boundary layers they 
have also been shown to modulate the near- wall turbulence, see for example [Hutchins~fcl 
Marusic (20076); Mathis et al. (2009), suggesting that they may play an important role 



in flow dynamics across a range of scales. 

In Couette flow, structures reminiscent of VLSMs have long been observed in the 
core through DNS of turbulent plane Couette flow (Lee fc Kim||1991 Bech et 0^1 11995). 
Although some studies raised the concern that the structures were numerical artifacts. 



recent DNS at higher resolution and with longer box sizes (Komminaho et al. 1996 



Tsukahara et a/. [2006 ) have confirmed the existence of long streamwise alternating high 



and low speed streaky structures at the centerline. In experiments, VLSMs were first 
identified through observations of a noticeable peak in the Fourier energy spectrum of 



the turbulence intensity at low frequencies (Komminaho et al. 1996 Kitoh fc Umeki 



2008). The Couette flow experiments of Tillmark fc Alfredsson (1998) found further 
evidence of very long structures in the form of long autocorrelations Ruu{^) and two- 
point correlations Ruui^x) as well as periodic variation of spanwise correlations Ruui^z) 
in the core. The streamwise extent of these correlations was longer than those generally 
seen in other wall-bounded flows. Komminaho et al. (1996) also found that in contrast 
to other flows, the streamwise correlations for Couette flow are larger at the center than 
near the wall. At channel center the zero cross distances of Ruuij) and Ruui^x) have 
been observed to be three times that of the corresponding structures in Poiseuille flow 
( Kitoh et fl/.||2005 ). This makes Couette flow an ideal candidate to test the applicability 
of a streamwise constant model. 

Streamwise constant, kx = 0^ perturbations to the LNS also produce the largest input- 
output response for both laminar (Bam ieh fc Dahleh 2001 Jovanovic fc Bamieh 2001 
Farrell fc Ioannou|l9936 1998} and turbulent ( del Alamo fc Jimenez 2006) base velocity 



profiles. In addition, streaks of streamwise velocity naturally arise from the set of initial 



conditions that produce the largest energy growth (Butler fc Farrell 1992; Farrell fc 



loannou 1993a), namely streamwise vortices. Even in linearly unstable flows, studies 
have shown that the amplitude of streamwise constant structures can exceed that of 
the linearly unstable modes ( Jovanovic fc Bamieh||20Q4 Gust avsson|| 1991 ). For channel 
flows, Bamieh fc Dahleh (2001) explicitly showed that streamwise constant disturbances 
produce energy growth on the order of Re^ whereas streamwise varying disturbances 
grow as a function of Re^ . 

In the present work we employ a streamwise constant model based on the previously 
discussed experimental and analytical evidence of the importance of streamwise homoge- 
nous features. This so-called two-dimensional, three (velocity) component, henceforth 
2D/3C, model for plane Couette flow is simulated under small amplitude Gaussian forc- 
ing. The results demonstrate the ability of this model to capture some important features 
of fully developed turbulent flow. In particular, it is demonstrated that: (1) the nonlinear 
terms in the 2D/3C model capture the momentum redistribution mechanism involved 
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Figure 1. Flow geometry. Streamwise and spanwise boundaries are periodic, the bottom wall 
is stationary and the top wall moves in the x-direction with a velocity Uw The channel half 
height is denoted 6 and the full channel height is denoted h. 

in creating the shape of the turbulent velocity profile, (2) a stochastically forced 2D/3C 
model can reproduce the appropriate turbulent mean velocity profile and Reynolds num- 
ber trends, and (3) this model produces amplification of small disturbances that is con- 
sistent with input-output studies of the LNS. The work is organized as follows: the next 
sections of the paper describe the model and simulation approach. Results and discus- 
sion follow, including a comparison between the model and a DNS dataset, before final 
conclusions. 

2. The 2D/3C Model 

The 2D/3C model discussed herein is obtained by setting streamwise (x-direction) 
velocity derivatives in the full NS equations describing Couette flow to zero ( |Bobba| 
2QQ4| ). This can be thought of as a projection of the NS into the streamwise constant 
space. One can explicitly show that for Couette flow this 2D/3C formulation also results 
in a system with zero streamwise pressure gradient. 

The velocity field is decomposed such that u = [U -\- u'^^; V + v^^^ W + ujsw]''> where 
U = U{y) = V = W = is the laminar Couette flow and {Ug^^Vg^^Wg^) are the 
corresponding time-dependent deviations from laminar in the streamwise constant sense. 

The flow geometry is shown in Figure [l] The Reynolds number employed is Rcu) = 
where is the velocity of the top plate, h is the channel height and u is the kinematic 
viscosity of the fluid. All distances and velocities are respectively normalized by h and 
Uw. In the sequel we will use {Ug^^v^^^Wg^) to denote (^^f^, ^jf^^i and explicitly 
indicate the scaling only in the figure labels. 

A stream function 

^SW Q^^ ^SW Qy 

forces the appropriate 2D continuity. This results in the following model: 

du',^ ^ djj du'^^ _ di^dl^ ^ dip du'^^ ^ I 

dt dz dy dz dy dy dz Re^ n 

d/^ip _ dip OAiP ^ dip dAiP ^ 1 * 

dt dz dy dy dz Re^ 

No slip boundary conditions at the wall and periodic boundary conditions in the spanwise 
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direction are applied (without loss of generality they can also be used for the stream 
function equation since v'^^ = w'^^ ~^^§f~^~^^V^~ const). This model 
retains many of the important flow features lost in a purely 2D model by maintaining 



all three velocity components. Equations (2.1) are an improvement over linear models 



because it is hypothesized that it is the nonlinearity in the u'^^{y^z^t) equation that 
provides the mathematical mechanism for the redistribution of the fluid momentum. This 
redistribution creates larger streamwise velocity gradients in the wall-normal direction 
and changes the plane Couette velocity profile from linear to its characteristic turbulent 
"S-shape". Meanwhile, the important features of the LNS are maintained. Linearization 



of (2.1) around the laminar profile produces a non-normal operator with a coupling term 



analogous to the one in the LNS. 

The laminar flow solution of Equation (2.1) was previously shown to be globally, that 
is nonlinearly, stable for all Reynolds numbers (Bobba et al. 2002), and therefore the 



laminar flow constitutes a unique solution. Consequently any transition mechanisms as- 
sociated with bifurcations, escape from the basin of attraction of the laminar solution or 
the like are not possible. So, any complications associated with these nonlinear phenom- 
ena can be eliminated from the analysis of these particular equations. Global asymptotic 
stability of the laminar solution also implies that without forcing, perturbations will 
eventually decay, in agreement with the results of 'Orlandi fc Jimenez] ( 1994[ ) who found 



that after an initial perturbation a 2D/?>C model decays (back to laminar) with time. 
The fact that one can analytically prove that the unforced 2D/3C model has a unique 
solution suggests that it is far more analytically tractable than NS. We do not pursue 
analytical studies of the 2D jZC model in the current work, but instead concern ourselves 
with showing the applicability of the model in describing important features of the flow 
field. However, the fact that global statements about these equations can be made implies 
that future analytical studies are promising. 

As with any model, there are assumptions built into the 2D jZC model, and it is im- 
portant to understand how these relate to the physical phenomena associated with tur- 
bulent flows. Most obviously, small scale, three-dimensional turbulent activity, including 
the specifics of several structures that are known to exist in the full flow, is not captured. 
While this makes appropriate scaling relationships more difficult to determine, it does 
not diminish the potential of the model for predicting and understanding key aspects of 
turbulence in plane Couette flow. The challenge lies in extending the 2D jZC model to 
incorporate aspects of the streamwise variation associated three-dimensional turbulent 
flow. 

2.1. Modeling Framework 

No model is a perfect representation of reality. In addition to modeling assumptions, 
parameter errors or external influences on the system in question are often ignored. 
Inaccurate parameter estimates or linearization of a nonlinear system may change the 
model's ability to predict behavior. Environmental conditions that affect (or disturb) the 
system may also play an important role in its dynamics. This role is not captured by a 
typical model. 

Robust control theory has historically been used to analyze models in the presence of 
such modeling errors (^uncertainty') [Poyle et al (1991); Zhou et al. (1996). One typically 



represents all of the uncertainties using an uncertainty operator A. The block diagram 
of Figure [2] is then used to depict a model subject to this uncertain set A. Generally 
robust control tools provide a bound on ||A||, below which a desired property can be 
maintained. Robust control tools do not require a detailed model of the particular un- 
certainty. This makes them appealing in situations where there are unknown (or hard to 
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A 









Figure 2. A robust control block diagram for a model 
subject to uncertainty. Generally a norm bound on A 
specifies the amount of uncertainty that a model can 
have before a desired property is lost (i.e., if the model 
is stable for ||A|| ^ 1, this implies robust stability). 



Zero Mean 
Noise 



2D/3C Model 



Figure 3. The approximation for il- 
lustrating the 2D/3C model's lack of 
robustness. The zero-mean noise as an 
approximation for the modeling errors 
and uncertainty. The noise acts as an 
additive "uncertainty" at each time 
step. 



model) environmental influences on the system, or when one can only specify the range 
on a parameter, rather than an exact value. However, since the uncertainty is generally 
specified through a bound that includes the worst case scenario, the results of this type 
of analysis may be very conservative. One way to mitigate this is to 'structure' or shape 
the uncertainty, a process which relies on some understanding of the implicit modeling 
errors. 

In the context of a system comprising a wall-bounded shear flow, many disturbances 
can be modeled through the A block in Figure |2j These sources of modeling errors (uncer- 
tainties) can arise from assumptions on the boundary conditions or unmodeled dynamics. 
External sources of model uncertainty that are not captured in the NS equations, include 
phenomena such as acoustic noise and thermal fluctuations in an experiment, or the 
build up of numerical error in simulations. In addition the uncertain set includes terms 
excluded by the modeling assumptions, namely the modes in the 2D/3C model, 
or the nonlinear terms for the LNS model. See Bobba (2004) for a full characterization 
of the types of uncertainties present in shear flow problems. Obviously the latter class 
of perturbations are strictly bound to satisfy the NS equations, while the former are less 
constrained. Distributed wall roughness (i.e., surface imperfections present in any real 
surface), wall vibration, imperfect alignment of the walls or other parameter estimates 
may be captured through either stochastic or other forcing. 

In the present work, the framework of robust control is employed in a nontraditional 
manner. Instead of providing an upper bound on ||A||, (i.e. a robustness guarantee) 
we describe the extent to which the laminar flow state is 'fragile' (i.e. unable to be 
maintained in the face of infinitesimal disturbances). One can think of this as an inverse 
robustness (or 'fragility') problem, i.e. a discussion of a lack of robustness. In order to 
study the disturbance response of the 2D/3C model the system of Figure [2] is abstracted 
into the simplified setting of Figure [s] We further simplify by linearizing the A?/^(?/, z^t) 
equation which is equivalent to recognizing that advection terms in the stream function 
equation play a lesser role in redistributing momentum. The forcing and henceforth 
iIj are constrained to be small such that the nonlinear terms are at least an order of 
magnitude smaller than the linear ones in all cases studied here. This small amplitude 
noise assumption is very important in the development of this work because of the focus 
on the effect of small amplitude disturbances on a fragile system and because larger 
amplitude forcing can change the dynamics of the model. 
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For all of the numerical studies described herein we simulate 
dz dy dz dy dy dz Re^ 



dt 



1 



(2.2) 



with the same boundary conditions as in Equation (2.1). A short simulation study com- 



paring low-order streamwise velocity statistics supports the use of linearization in the ijj 
equation. 

Approximating the full 31) system using the interconnection in Figure [2] would involve 
nonlinear mixing modes. In order to approximate the range of frequencies associated 
with the full 3D system in the framework of Figure [3j zero mean stochastic forcing was 
applied to the 2D jZC model. In particular, the inputs du{y,z^t) and d^ijj^z^t) in (2.2) 
are small amplitude and Gaussian, as in Gayme et a/. | ( [2009 ). The input amplitudes are 
defined using the standard deviation, (Jnoise- Note that under these assumptions there is 
no coupling from the streamwise components {u) back to the cross-stream components 
(the All) equation). The plausibility of modeling the type of disturbances common to 
experimental conditions in this manner is confirmed by results from stochastic forcing of 
the LNS equations, which leads to fiows dominated by streamwise elongated streaks and 



vortices that are strikingly similar to those observed in experiments (Farrell & loannou 



19936), as well as by the results of the simulation study discussed in Section 4.3 Further 
development of the model would be required to address this effective feedback mechanism. 



3. Approach 

Time-dependent simulations of the full coupled system \2.2\ were carried out using a 
basic second-order central difference scheme in both the spanwise {z) and wall-normal {y) 
directions. Periodic boundary conditions in z and no-slip boundary conditions in y were 
applied. Simulations using the spectral methods of Weideman & Reddy (2000) were also 
performed for comparison. The pseudospectral simulations employ a Chebyshev inter- 
polant for the wall- normal direction and a Fourier method for the spanwise derivatives. 
The aspect ratio in all of the simulations was greater than 12 to 1 (spanwise to wall- 
normal) in order to eliminate box size effects; specifically the usual computational box 
size was Ly x = h x 12.Sh with 75 x 100 grid points. A spanwise extent of 12.Sh was 



selected to provide a direct comparison to the full field DNS data from Tsukahara et al. 
( 2QQ6| . 

In this study, the response of the streamwise velocity, li'^^, to forcing of the cross- 
stream velocity components, v^^ and w^^, was examined. A forcing input of zero mean 
small amplitude Gaussian noise evenly applied at each y-z plane grid point was selected 
for d^p. The other input forcing, d^, was set to zero based on previous studies of the LNS, 
which showed that the response to streamwise body forcing is significantly smaller than 
the response to spanwise/wall- normal plane forcing ( Jovanovic & Bamieh 2005[ ). These 
studies used an order of magnitude argument to conclude that the difference in response 
scales as Furthermore, it is energy redistribution by streamwise vorticity (i.e. A^) 
that is thought to be the primary effect governing the shape of the turbulent velocity 
profile ( [Hamilton et a/.|1995D . The response in the streamwise velocity component to this 
forcing may have a nonzero mean because of the nonlinearity in the u equation. 

The different discretization techniques naturally provide a comparison of different noise 
forcing distributions. For example, the Chebyshev grid results in a higher concentration 
of noise forcing near the walls. Throughout the present work it is assumed that significant 
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numerical errors are not introduced by the methods of discretization, i.e. the introduction 
of significant noise arises only through the d terms of Equation (2.1). 

The time evolution of Atp in Equation (2.2) can be seen to be a stochastically forced 
heat equation, i.e. a linear stochastic partial differential equation which can be solved 
analytically (see for example Swanson^ (2007) or Luo^ (2006) and the references therein). 
This is not pursued here because a simulation is a much simpler way to demonstrate 
the efficacy of the model. An exposition on Ito calculus and Wiener chaos expansions is 
beyond the scope of this paper. Future work may involve pursuing analytical solutions 
to both the linear approximation to and the full nonlinear system (2.1). 



4. Results and Discussion 

The results will be divided into three main sections. We begin by analyzing the DNS 



data of Tsukahara et al. (2006) in the light of the 2D/3C model and confirming the 
extent to which the assumptions of the 2D/3C model can be adduced through this 



data. Following that, a time independent version of Equation (2.1) is studied to verify 



the implicit model filter between ip and u. Finally, results from full simulations of the 



time-dependent Equations (2.2) are presented and compared to the DNS data. 



4.1. Comparison of DNS data with 2D/3C Modeling Assumptions 



Full details of the DNS dataset can be found in Tsukahara et al. (2006); a brief review 



of key aspects is given here. Three Reynolds numbers were considered, Re^ = 3000, 8600 
and 12800, all with computational domain size x Ly x = 44.8/i x h x 12.8/i, 
1024 X 96 X 512 grid points, and a sampling time (^^) of 91. The fourth-order finite 
difference scheme proposed in Morinishi (19951) was employed for the x and z directions. 



A second-order finite difference method was used for the y direction. 

The friction coefficient, Cf = 9.59 x 10~^, is somewhat higher than in other studies. 



such as [Robertson fc Johnson (1970). Filling this friction factor into the relationship 
developed by Robertson]! 1959 ) , 



G 



where C / 



(4.1) 



with Tw used to denote shear stress at the wall, leads to an experimental constant G = 
0.199. Other values reported in the literature include G = 0.19 and G = 0.174 both 



from Robertson ( 1959 ) based on the data of Reichardt and Robertson respectively and 
G = 0.182 from the experimental study of |E1 Telbany fc Reynolds] (|1982D. 



The turbulent mean velocity profiles, turbulence intensities, Reynolds stresses and 
budgets of u^u^^ from this DNS show good agreement with the experimental results of 



Tillmark (1995) and the spectral DNS study of [Komminaho et al. (1996) which used a 



larger box. The two-point correlations in u indicate that the box lengths used in both 
the streamwise, Ruu{^x)^ and spanwise, i?nn(Az), directions are sufficient to eliminate 
any boundary condition-related spurious effects. 

In what follows, a streamwise constant projection of the DNS data is approximated 
through a streamwise {x) average over the box length, which will highlight streamwise 
coherence of the order of the box length. The x-averaged DNS data is denoted u^^^ = 
('^Xave ^ ^(y)^^Xave^'^Xave) distiuguish It from truc streamwise constant data. Time 
averages are indicated by an overbar, (•). 

The ratio of the energy contained in the x-averaged DNS to that of the full field provides 
a quantitative measure of the extent to which the DNS data can be approximated as 
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Component 


Total Energy Norm 

Ml 


Percent of Total Energy 
in x-averaged Norm 


u 


0.5334 


99.1 


u-U 


0.1686 


90.2 


V 


0.0279 


19.0 


w 


0.0412 


15.0 



Table 1. Energy content in the x-averaged DNS velocity components at Rcw = 3000. 



streamwise constant. For this comparison the squared 2-norm is used to approximate the 
energy in each 2-dimensional x-averaged velocity component 

m?= I C P{y,zfdydz 
JZ Jo 

where Z is the spanwise extent, = 2^2 — ^1 is the spanwise distance between z grid 
points and trapezoidal approximations are used for the inhomogeneous y grid. 

Table [l] shows the total energy (based on the full DNS field at Re^ = 3000) and 
the percentage contained in each of the x-averaged velocity components (li, w) as 
well as in the deviation from laminar (denoted u — U). This latter quantity is most 
representative of the energy associated with the differences in the mean velocity profile 
for a turbulent versus a laminar flow. The computations show that x-averaged streamwise 
velocity contains 99% of the {u) energy, whereas the corresponding deviation from laminar 
contains 90%. As expected, the x-averaging results in a larger loss of information in the 
spanwise and wall-normal velocity components. 

An examination of the DNS streamwise velocity field at =29, close to the outer edge 
of the region affected by the near-wall cycle, reveals the signature of streamwise elongated, 
large scale streaks in the streamwise/ wall- normal plane of the full field (Figure [4|. These 

and [5(b)] which depict contour plots of the deviation 
from laminar flow, u'^ = Ux^^^ — U ^ when averaged over 25% of the streamwise field 
and the full field respectively. Clearly, increasing the averaging length acts as a filter on 
structures of different streamwise extent. The average over the full box length retains 
strong evidence of structures across the entire spanwise/wall- normal plane. In particular, 
the strongest signature near the wall is in qualitative agreement with the near-wall model 
of energetic structures centered around y+ ^ 15 with a statistical diameter of y^ ~ 30. 
Another important feature of Figure |5(b)| is that the maximum deviations from laminar 
flow, which are out of spatial phase with one another, top to bottom, are associated with 
large-scale rolling motions which reach across the channel height. 

The above analysis shows that there is good agreement between the DNS data and our 



streaks are also visible in Figures 5(a) 
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Figure 4. An z-x plane contour plot of the streamwise velocity, u, from the DNS field, (bottom 
up view) at y~^ = 29. Light colored contours denote regions of higher velocity and dark contours 
indicate lower velocity regions. 




(a) Average over 25% of streamwise box (b) Average over full streamwise box 



Figure 5. y-z plane contour plots of the x-averaged (as an approximation for streamwise 
constant) DNS deviations from laminar (ux^^^) (a) averaged over 25% of the streamwise box 
length and (b) averaged over the full streamwise box length. 



assumptions. In the next section this data is used to suggest a time-independent model 



for ip{y^z) in order to study a steady state version of the streamwise velocity in (2.1). 
This is followed by simulation of the full system (|2.2[). 



4.2. Time Independent Ux{y^ z^t) Equation 



The response of Equation 2.1 to a stream function that is independent of time, ^l^ssiv-, z)^ 
is of interest for two reasons: (1) the forced solution for the streamwise velocity permits 
investigation of whether the 2D jZC model filters an appropriately shaped ipss{y^z) to- 
wards the expected shape of the turbulent velocity profile, and (2) it gives insight into the 
mathematical mechanisms that create the momentum (energy) transfer that generates 
the blunted profile. The analysis constitutes a weakly nonlinear analysis, in which the 

time- independent forcing takes the form iIjss = i^sso + ^i^ssi + 

In [Barkley fc Tucker man| ( |2QQ7j it was shown that laminar-turbulent flow patterns in 
plane Couette flow could be reproduced using a stream function of the form ?/^(y,z) = 
ipo{y) -\-ilji{y) cos{kzz) +'02(^) sm{kzz). We use this study as guidance but set the zeroth- 
order term t/^q to zero because a nonzero i/jq produces a nonzero- mean spanwise flow w^^ , 
which is not representative of the velocity field we are interested in studying. The DNS 



field (Tsukahara et al 2006) was also used as a guide to ensure that the first-order 



term ipssi as well as the corresponding wall- normal and spanwise velocities, respectively 
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and 



contained representative features. For ease of computation and analysis 



a simple analytic model for ij^ssi [y-, z) was selected, 
which matches the boundary conditions 



namely a doubly-harmonic model 



This corresponds to 
27r 



27r 



esin"^ (Try) cos ( -^z 



(4.3) 



27r 



The size of the perturbation, e, is a free variable to be explored, while the spanwise 
wavelength, A^, is fixed to a value determined using the DNS data. 

Streamwise averages of both v{x^y^z) and w{x^y,z) from the DNS data permit an 
estimate of ^l^ss{y^ (to within some constant), for that particular field. A contour plot 
of the approximation based on w'^^^^ z) is shown in Figure 6(a) The value of 



l.Sh 

was chosen to match the results from a Fast Fourier Transform (FFT) of this data while 
maintaining the same box size (12.8/i) employed for the DNS. This value is also in the 
range of the spanwise wave number corresponding to maximum amplification of the 
linear operator (optimal spanwise spacing), kz G [2.8,4] (A^ G [1.6, 2.2]/i), reported in 
the literature ( Farrell fc Ioannou|[T9936{ [Butler fc Farren||1992[ |Gustavsson]|1991| ). An 
initial perturbation amplitude of s = 0.00675 was selected based on the approximate 
values obtained by integrating v'^^^{y^ z) and w'^^^{y^ z). The estimated amplitude is very 
small, in agreement with the idea of using a nominal model plus an uncertainty which 
is amplified through the coupling in the linear operator in a manner that is described 
and quantified in studies such as Trefethen et al. ( 1993') and ' Jovanovic fc Bamieh| ( |2005 ) . 

A contour plot refiecting these parameter values is provided in Figure [6(b)[ It shows 
good qualitative agreement, in particular with the region of strongest signal in the DNS 
streamwise average (Figure [6(a) ). The latter wall-normal variation is complicated (and 
Reynolds number-dependent), but a simple harmonic variation gives a reasonable rep- 
resentation. The velocity vector field implied by Equation (4.3) is consistent with low 
speed fiuid being lifted up from the stationary wall and higher speed fiuid being pushed 
down from the moving wall, and as such supports the notion that the mechanisms of 
interest can be modeled using a single harmonic in both y and z. 

The stream function of Equation |4.3| with the selected e and A^ was applied to the 
time- independent form of equation |2.2[ yielding 



d^ss d ^ dips 
dz dy dy 



1 



dz Re,. 



d^ss dU 
dz dy ' 



(4.4) 



A contour plot of the resulting u^^^ (?/, z) is depicted in Figure [7[ This fi ffure shows a 
'^swss (^' ^) with near- wall rolls that are out of spanwise phase with one another similar to 
those seen in the x-averaged DNS data of Figure [5'(b)| The increased coherence associated 
with the steady state model relative to the DNS data manifests as an increased variation 
in the deviation from laminar (amplitude of the surface) particularly at the center of the 
channel. This effect is emphasized through comparison of the surface plots of Figures 
a) and |8(b)] Note the different vertical axis scales for the two plots. 
Averages across the span of u^^ (?/, z) for e = 0.00675 as well as for four additional £ 



values are compared to a similar average of u'^ave the DNS in Figure 9(a) Clearly 

using i/jss from (4.3) as an input to (|4.4[ ) produces streamwise velocity profiles whose 
shapes are consistent with u^^ ^ u — U from the DNS. The peaks are, however located 
at different wall-normal positions. An amplitude that exactly matched both the magni- 
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2 4 6 8 10 12 2 4 6 8 10 12 

1 ^ 



(a) ipxave Estimated from w^^^^ (b) ipssiy, z) = 0.0675 sin^ {ny) cos 

Figure 6. (a) Contour plot of the x- averaged DNS (as a streamwise constant, 
2D/?>C approximation) spanwise velocity deviations integrated to obtain the stream func- 
tion, 'ipxaveiy^^) — ~~~§§^- (^) Contour plot of the simple harmonic model for 
^l^ssiy, z) — 0.00675 sin^ (vry) cos {y^z) with amplitude and wavelengths that approximate DNS 
data. 




2 4 6 ^ 8 10 12 2 4 6 8 10 12 



(a) DNS Data (b) u'sw^^ predicted from ij^ss 

Figure 7. Contour plots of (a) x- averaged DNS data and (b) the 2D/3C (stream- 
wise constant) velocity deviations, u'avess^ obtained using steady state estimate 
i/jss{y, z) — 0.00675 sin^ (tt^) cos ( ^ 3143 ^) ; both plotted with the same contour levels. 



tude and location of the DNS peaks was not found even when different values of kz were 
studied. This is not unexpected because of the simplicity of the wall-normal variation in 
the steady state model, as well as the streamwise constant and steady state assumptions 
(clearly the full turbulent field is neither streamwise constant nor time-independent). 
However, this type of model clarifies the nonlinear role of cross-stream flow features in 
redistributing energy in the flow field. These results suggest that the phenomenon that is 
responsible for blunting of the velocity profile in the mean sense is a direct result of the 
interaction between rolling motions caused by the y-z stream function and the laminar 
profile. In other words, this study provides strong evidence that the nonlinearity needed 
to generate the turbulent velocity profile comes from the nonlinear terms that are present 
in the u'g^{y^z^t) equation of the 2D/?>C model (2.1). 

The magnitude of forcing applied to the system is reflected in the amplitude of ^pssiv^ ^) 
which in turn affects the friction Reynolds number, Rcr through the skin friction arising 
due to the resultant mean velocity gradient at the wall. Increasing the amplitude (s) in 
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(a) u'x^^^ DNS data (b) u'sw^^ estimated using ipssiy^z) 

Figure 8. Surface plot of velocity deviations (a) u^^^^ {y, z) from DNS Data at Rcw = 3000 and 
(b) u'swssiv^ obtained using steady state estimate il^ss{y,z) — 0.00675 sin^ (tt^) cos ( ^ §143 z) 
at Rcw — 3000, note the z scale difference between (a) and (b) 
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(a) Amplitude Variation 




Amplitude 

(b) Velocity Gradient at the Wall 



Figure 9. (a) Variation of the 2Z)/3C (streamwise constant) velocity deviations, u'^^^^ with 
A] estimates are obtained using steady state estimate ip{y,z) = A sin^ (iry) cos ( g^^g z) (b) 



Variation in the velocity gradient at the wall 



dy 



with ipss amplitude, A. 



Equation 4.3 is analogous to increasing the magnitude of the model uncertainty. These 
effects are emphasized in Figure |9(b)| which provides a plot of e versus the velocity 



gradient at the wall. This behaviour will be discussed in more detail in Section [43 

4.3. Time Dependent 2D/3C Model 



Time-dependent simulations of (2.2) were carried out using the basic second-order central 
difference scheme and pseudo- spectral approaches described earlier. Table [2] lists the 
Reynolds number and forcing amplitude combinations considered. The window used for 
time averaging was At = 100000^. 

The initial simulation (Case 1, in Table [2]) was carried out at Rcu, = = 3000 with 
d^{x, t) drawn from a zero mean Gaussian distribution with standard deviation (noise 
amplitude) a noise = 0.01 applied at every point in the mesh. 

A comparison of Figures |5(b)| and |10(a)| shows that contours of constant streamwise 
velocity deviation from laminar from the DNS and the 2DjZC simulation are in good 
qualitative agreement. In particular, the spanwise offset in spatial phase between peaks 



14 D. F. Gayme, B. J. McKeon, A. 



Papachristodoulou, B. Bamieh & J. C. Doyle 




(a) 




Figure 10. (a) Contour plot of t)obtained from the 2D/3C model for Case 1 with 

the same contour levels as in Figure l7| (b) The surface plot corresponding to (a). 



Table 2. Computation Details 



Case Reynolds Number 


O'noise 


Ly 


X Lz 


Ny X Nz 


Squared Norm of 
the Noise Input 


1 


3000 


0.01 


h X 12.8/1 


75 X 100 


0.0565 


2 


3000 


0.0125 


h X 12.8/1 


75 X 100 


0.0882 


3 


3000 


0.004 


h X 12.8/1 


75 X 100 


0.009 


4 


8600 


0.004 


h X 12.8/1 


75 X 130 


0.0092 


5 


12800 


0.004 


hx ^ 


- 16.5/1 


75 X 130 


0.0092 


6 


12800 


0.001 


hx ^ 


- 16.5/1 


75 X 130 


5.77e - 04 


Spec 1 


3000 


0.001 


hx ^ 


- 14.5/1 


40 X 81 




Spec 2 


3000 


0.002 


hx ^ 


- 14.5/1 


40 X 81 





from top to bottom is reproduced. While the dominant wavelength from the 2D/3C 
simulation is somewhat longer than the ^ 1.8 of the DNS (frequency analysis of u'^^ 
indicates that most of the energy resides in wave lengths between 4 ^ ^ 6.1) there is 
also a significant contribution from A^ ~ 2. There is noticeably better agreement between 
the time-dependent model and the DNS (Figure p-0 (b) | and 8(a) ) than for the steady state 
analysis (Figure [8 (b)| ), likely a consequence of the broadband stochastic, i.e. less coherent 
and time-dependent, forcing. 

4.3.1. Mean Velocity Profile 

Averaging the streamwise velocity field obtained for Case 1 in Table|2]leads to the mean 
velocity profile (i.e. iT^) shown in Figure 11(a) The mean profile can also be plotted in 



A Streamwise Constant Model of Turbulence in Plane Couette Flow 



15 























































f 


















J 














































f 


















J 


























2D /3C Model, Oa 
•--DNIS data Rt„ = 
Laminar Profile 


















3000 















0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



(a) Mean Velocity Profiles 



-DNIS data Rs^ = Si. 939 
-2D/3(7Mbdel, Gase 2 
-2£)/3(7 Model, OaEe 1 



(b) Mean Velocity Inner Units 



Figure 11. (a) Comparison of mean velocity profile from the 2D /?>C model Case 1, Usw{y, z, t), 
with u{x,y, z,t) from DNS. (b) Inner scaled velocity profiles comparison of Cases 1 and 2 to the 
DNS data with Rer ~ 52 for all data sets. 



inner units (Figure 11(b) ) with the use of Equation 4.1 (with G = 0.1991 from Tsukahara 
et al. ( 2QQ6| )) to estimate the friction velocity, u^- There is good agreement between the 



DNS and the Case 1 simulation, even with the assumption of a friction velocity that 
corresponds to the full flow. However, it is clear that below ~ 20 the 2D jZC model 
underestimates the expected velocity profile (maximum error 7.4%), and above that it 
overshoots it (maximum error 2.4%). There are two obvious first-order interpretations of 
these discrepancies. First, for cases 1-6, the noise is modeled as being evenly distributed 
across the grid while in reality the noise is likely higher in the buffer region due to the 
proximity of the wall, and lower in the overlap layer. An improved noise model might 
improve the agreement. A second interpretation is that further from the wall the flow is 
better modeled by the streamwise constant approximation, while streamwise variation is 
more important in the dynamics of the near-wall region (in agreement with the known 
variation of the spectral distribution of streamwise energy in the full flow). 

A second (constant) noise amplitude at the same Reynolds number. Case 2, is also 



shown in Figure 11(b) The agreement with the DNS is certainly improved below y'^ = 20 
(maximum error 6.19% at y~^ = 19), but at the expense of larger error further from the 
wall (~ 5 — 6% between 20 < y^ < 30). These results further support the idea that 
a non-uniform noise forcing with increased noise near the wall versus that at channel 
center may more accurately reflect the conditions in a real flow field. This idea is further 



explored in Section [4. 3. 3[ 



It should be noted that Ur can also be computed directly from the velocity gradient at 
the wall. In both cases Rcr is underestimated by around 10% compared to the estimate 
from ( |4.1| ). Because of the limited number of points near the wall, the friction relationship 
from the full flow was preferred, with the understanding that this would only be correct 
if the 2D/3C model with a^oise exactly reproduced the mean flow behavior. 

4.3.2. Reynolds Number and Noise Amplitude Trends 

Four additional Reynolds number and (Jnoise amplitudes pairs were considered. The 
details for each of the cases 3-6, along with the computational domain and spatial res- 
olution, are provided in Table |2] Respective values of the norm || • |p, as computed in 



Equation (4.2), of the noise input computed over the box are also reported, since this is 



a more appropriate measure of the forcing when the box size varies. 
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It is useful to introduce a normalized version of Equation (2.1) through the chan ge o f 
variables t = and ^ = Re ip. This creates new expressions for the forcing in (2.2), 
Vu = Redu{= 0) and = Re'^d^. The expression = Re'^d^ indicates that an 
increase in noise produces a similar effect to an increase in Reynolds number (actually 
^/Re) , as observed in the increased deviation from laminar observed with increasing noise 
amplitude in |9(a)[ This is especially clear when considering the variation of Rer because 
an increase in noise amplitude directly corresponds to increased velocity gradients at 
the wall due to the no-slip boundary conditions. Increased profile "blunting" with both 
increasing anoise (noise input energy) and Reynolds number in cases 2-5 can be observed 
in figure 12(a)[ 



For the higher Reynolds numbers (but constant noise amplitude) in Case 4 and 5, 
Figure |12(b)| sho ws a worsening agreement in inner units with the DNS data from 
Tsukahara et al ( 2QQ6[ ) at Re^- = 128.5 and Re^- = 181.3, respectively. The underes- 



timation below ^ 30 (in the buffer layer) is more pronounced, but the agreement 
above y~^ > 30 remains of similar magnitude (max error ~ 4.94% for Re^- = 128.5 and 
8.39% for Rer = 181.3). We hypothesize that this worsened agreement may be repre- 
sentative of the increasing scale separation with increased Reynolds number. Near- wall 
motions that can effectively be modeled as streamwise constant at low Reynolds numbers 
have an increasingly short streamwise wavelength relative to the motions that scale with 
outer length scale S. That the zero-error location consistently occurs around y~^ = 20- 
30, commonly thought to be the upper boundary of the buffer layer, is consistent with 
this scale separation argument. For the same reason, the lack of model resolution in the 
near- wall region will be exacerbated with increasing i?e^. In robust control terms, this 
points once again to an increase in the model uncertainty near the wall versus the channel 
center. Once again, a better uncertainty model could be accomplished through the use 
of a 'structured uncertainty' which would include an increase in anoise in the near-wall 
region. 

4.3.3. Varying Noise Amplitude and Distribution 

A preliminary effort to introduce a non-uniform distribution of noise was carried out 
by repeating the simulation using a pseudospectral scheme with a Chebyshev interpolant 
for the wall-normal direction. This scheme naturally produces increased noise near the 
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Figure 13. (a) Comparison of u'^ versus from 2D/?>C Model using Chebyshev spacing in y 
with DNS data at Rcw — 3000 based onG — 0.1991 [Rct ~= 52) (b) Comparison of Usw{y, t) 
from 2D l?)C Model at for Case 1 [Rcw — 3000 with Gnoise — 0.01, energy — 0.0882), and Case 
6 (Rew — 12800 with Gnois& — 0.001, energy — 5.77e — 04) same grid and box size. 



walls. Cases Spec 1 and Spec 2 in Table [2] are two such simulations, both at Re^j = 3000, 



with (Jnoise = 0.001 and (Jnoise = 0.002 respectively. Figure 13(a) shows the resulting 
mean velocity profiles. Clearly the noise level is too low for Spec 1, however for Spec 2 
the maximum error occurs in the buffer layer and is of the order 5-6%. The results of 
the spectral simulations indicate that by further noise shaping one could improve the 
agreement throughout the profile and across a range of Reynolds numbers. 

As previously discussed there is a strong relationship between the friction Reynolds 
number and anoise- As an illustration of this. Figure [T3(b)| shows that one can obtain 
similar mean velocity profiles at two different Reynolds numbers simply by adjusting 
the noise amplitude, i.e. a higher Reynolds number requires a smaller (uniform) noise 
amplitude to develop a mean velocity profile that is similar to that of a lower Reynolds 
number case with higher noise amplitude. This result is consistent with observations 
of higher transitional Reynolds number associated with "quiet" experiments compared 
to ones with high background disturbance levels. Alternatively, a fixed amplitude noise 
produces a larger response (more blunting) at higher Re than lower ones because dis- 
turbance amplification increases with increasing Reynolds number. This example makes 
it clear that the noise amplitude and the friction Reynolds number are tightly coupled, 
while giving further evidence that Reynolds number dependent wall-normal shaping of 
the noise would be required to get a better model representation of the turbulent mean 
velocity profiles. 

4.3.4. Characterization of the (small) disturbance amplification 

The results described herein indicate that a very small amount of stochastic noise 
forcing limited to the cross-stream components produces a very large response which 
corresponds to a behaviour that is not a solution of the unforced equations. The ability 



of our model (Equation 2.1 ), which has a unique solution in the unforced case, to produce 
a new flow condition due to such a forcing supports the notion that the model is not 
robust to small disturbances/uncertainty. The potential for disturbance amplification is 
not new, in fact it comes directly from the features of the LNS previously discussed, 
however the creation and maintenance of the new flow state is different and cannot 
come through the use of a linear model. A simple characterization of the amplification 
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maintained through the forced response of (2.2), or the lack of robustness ('fragihty') of 
the system, can be formulated as follows. 



Defining the squared 2-norm of the streamwise component of (2.2) (i.e. H'^^g^yH ) to be 



the increase in streamwise energy from the base (laminar) fiow, a so-called amplification 
factor is given by, 

r« = (4.5) 

W^noise \ \ 

which is a measure of the output energy for a given input (noise forcing amplitude). 
is a nonlinear analog of the 'ensemble energy density' described in previous studies of 
the input-output response of the LNS, e.g. (Bamieh fc Dahleh|2QQl Jovanovic & Bamieh 
2005). Those investigations showed that the coupling between the Orr-Sommerfeld and 



Squire modes enables very large, Reynolds number dependent disturbance amplification. 
The amplification factor for case 3-5, which all have approximately the same input energy, 
are respectively Tu ^ 680, ^ 2200, and F^ ^ 2920. These trends are consistent with 
the low Reynolds number scaling trends based on the OSS equations. This agreement 
reflects the effective restriction of the streamwise constant assumption to amplification 
of the kx = modes (most amplified in the linear equation) as well as the source of the 
amplification in the 2D/3C model: a coupling in a similar linear operator, given by 



d 


f 




dt 







dU d 
dz 







(4.6) 



In this way computing F^^ from the simulation of (2.2 ) is analogous to studying the steady 
state nonlinear response to the most amplified 31) mode, (i.e. the kx = ^ mode). 

Equivalent amplification relationships between the cross-stream velocity components 
and (Jnoise could similarly be investigated. 

Structures with long streamwise coherence have long been shown to have a significant 
role in both transition and fully developed turbulent fiows. Based on these observations 
we study a streamwise constant projection of the Navier-Stokes equations for plane Cou- 
ette fiow, the 2D jZC model. Simulation of this model under small amplitude Gaussian 
forcing captures the turbulent mean velocity profile at low Reynolds numbers. Appropri- 
ate Reynolds number trends are also reproduced. 

A weakly nonlinear, steady state analysis demonstrates that 2D stream functions can 
produce appropriate mean velocity distributions when they are nonlinearly coupled to 
the 2D jZC streamwise velocity. This indicates that 'swirling motions' in the y-z plane 
produce features consistent with the mean characteristics of fully developed turbulence. 
It also provides evidence that the nonlinear coupling in the 2D jZC model is responsible 
for creating the well known characteristic 'S' shaped turbulent velocity profile. 

The use of small amplitude stochastic forcing as an input to the 2D jZC (nominal) 
model is based on ideas from robust control. Experimental observations are used to 
simplify the NS equations to form this nominal model and the noise forcing is used to 
capture both uncertain parameter values and unmodeled effects. The resulting forced 
2D/3C model allows one to isolate phenomena that can not be decoupled from a full 
simulation of NS while maintaining a sufficiently rich description of the physics that 
govern turbulent ffow. Our physics-based model should provide greater insight into the 
dynamics of the system than an empirical technique. Such a model may also allow better 
control design. 



The linearized 2D/?>C model (4.6) maintains the properties responsible for large dis- 
turbance amplification which have also been linked to subcritical transition. Maintenance 
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of these linear mechanisms is critical to the success of this approach. It is the combina- 
tion of these linear processes along with the momentum transfer from the two nonlinear 
terms in the streamwise velocity equation that enable the model to capture the blunted 
turbulent velocity profile. This line of inquiry provides a complementary perspective to 
transient growth and structurally based models, in that the 2D/3C model offers some 
improvement in analytic tractability at the expense of streamwise detail. The results are 
especially promising because the computational and analytical tractability of this model 
makes it well suited to higher Reynolds number studies. 

A natural extension of the present work would be the development of a more appropri- 
ate model for the noise distribution. It is common in the controls literature for a system 
to have a so-called structured uncertainty which is based on the physics of a particular 
system. In this work the limitation of noise to only the At/^ equation represents a first 
level of such an approach. Knowledge of the physics, for example that the near- wall re- 
gion is under- resolved in the 2D/3C model, is a first step. Numerical or experimental 
studies aimed at characterizing true spatial noise forcing patterns would further help in 
determining the correct model for noise distribution. 
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